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Abstract 


The  Cambridge  Atmospheric  Density  Numerical  Integration 
Program  (CADNIP)  is  a  completely  automatic  computer  program 
capable  of  determining  atmospheric  densities  from  an  analysis  of 
satellite  observations.  The  adopted  approach  consists  of  a  numerical 
integration  procedure  combined  with  a  differential  correction  scheme 
where  discrepancies  between  computed  and  observed  satellite  position 
and  velocity  are  reconciled  by  adjusting  the  assumed  atmospheric  model, 
thereby  yielding  corrected  or  refined  density  data. 

This  report  documents  the  latest  version  of  the  program  which 
includes  significant  refinements  and  improvements  incorporated  into 
CADNIP  under  this  contract. 
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Section  1 
Introd  action 

The  IBM  Cambridge  Advanced  SystemB  Department  has  been  actively 
engaged  in  the  design  and  development  of  automatic  computer  programs  ca¬ 
pable  of  determining  atmospheric  densities  from  an  analysis  of  satellite 
observations.  This  interest  (which  has  been  continuous  over  the  past  ten 
years)  has  resulted  in  two  totally  different  computational  approaches  to  the 
problem.  For  relatively  high  satellites,  where  the  effects  of  atmospheric 
drag  are  small,  a  successful  method* developed  by  IBM  consists  of  analyzing 
the  orbital  elements  of  the  satellite  over  long  periods  of  time.  In  particular, 
a  precise  knowledge  of  the  rate  of  change  of  the  period  may  be  used  for  den¬ 
sity  determination. 

For  low  satellites,  where  the  effects  of  atmospheric  drag  are  consid¬ 
erable,  an  alternate  approach  has  been  developed  by  IBM.  This  method 
consists  of  a  numerical  integration  scheme  combined  with  a  differential 
correction  procedure  where  the  discrepancies  between  computed  and  ob¬ 
served  satellite  position  and  velocity  are  reconciled  by  adjusting  the  as¬ 
sumed  atmospheric  model,  thereby  yielding  corrected  or  refined  density 

data.  An  important  advantage  of  this  technique  is  its  applicability  to 

2 

satellites  entering  the  decay  stage.  This  method  has  also  been  success¬ 
fully  incorporated  by  IBM  into  completely  automatic  computer  pr  ograms 
capable  of  computing  atmospheric  densities  from  an  analysis  of  optical 
or  electronic  observations  of  a  satellite. 

The  research  effort  reported  upon  herein  consisted  primarily  of  an 
attempt  to  improve  the  techniques  developed  under  the  latter  method.  These 
refinements  were  aimed  at  increasing  the  precision  of  the  computed  results, 
decreasing  the  required  computer  time,  and  incorporating  additional  flexi¬ 
bility  and  utility  into  the  program. 

As  a  result,  an  improved  version  of  the  Cambridge  Atmospheric  Den- 


aity  Numerical  Integration  Program  (CADNIP)  was  obtained  and  is  documented 
in  this  jTinal  Report.  Special  attention  should  be  given  to  Sections  5,  and  6 
along  with  Appendix  B  which  emphasize  the  aspects  of  the  research  effort 
performed  under  this  contract. 
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Section  2 

Procedures  Summar 


This  section  is  devoted  to  a  brief  description  of  the  procedures  used 
in  the  program. 

Given  a  set  of  satellite  observations  covering  a  minimum  time  span, 
it  is  desired  to  deduce  from  these  observations  the  density  of  the  atmosphere 
in  the  region  of  space  covered  by  the  satellite.  The  procedure  requires  an 
initial  estimate  of  the  vehicle's  orbital  parameters  along  with  an  estimate 
of  the  nighttime  exospheric  temperature  at  some  time. 

The  given  orbital  elements  are  converted  to  position  and  velocity,^ 

The  numerical  integration  scheme  is  then  called  for  tabulating  computed 
position  and  velocity  at  each  of  the  observation  times.  A  standard  differ¬ 
ential  correction  technique  is  then  used  to  reconcile  discrepancies  between 
computed  and  observed  position  and  velocity.  The  differential  correction 
scheme  actually  corrects  seven  quantities;  the  six  elements  of  position 
and  velocity  and  the  assumed  atmospheric  model.  Finally,  the  corrected 

position  and  velocity  are  converted  back  to  orbital  elements  for  output 
4 

pu  rposes. 

An  additional  feature  of  the  program  consists  of  a  pre.iminary 
adjustment  routine  which  attempts  to  improve  initial  estimates  of  certain 
parameters  most  likely  to  be  in  error.  This  feature  effectively  supplies 
the  differential  correction  routine  with  more  accurate  starting  conditions 
which  tend  to  increase  the  probability  of  convergence  and  also  to  speed 
up  the  procedure. 

Upon  completing  the  processing  of  one  set  of  observations,  the 
program  proceeds  to  the  next  set  of  observations,  with  the  results  of  the 
previous  set  being  used  as  the  initial  conditions  for  the  current  epoch. 


Section  3 


Mathematical  Considerations 

This  section'summarizes  briefly  the  major  mathematical  tech¬ 
niques  incorporated  into  the  program. 


3.  1  Numerical  Integration 

The  perturbative  forces  considered  for  this  study  are  the  earth's 
gravitational  field  and  air  drag. 


Geopotential  Representation 

The  earth's  geopotential  is  represented  by 


u 


.  m  .  m  — m 

(  C  cos  m\  +  S  sin  rru  )  P 
n  n  n 


where: 


r  a  radial  distance  of  the  satellite  from  the  center  of  the  earth 
X  =  subsatellite  longitude 

s  subsatellite  geocentric  latitude 

Cm,  Sm  are  the  earth's  zonal,  tesseral,  and  sectorial  harmonics 
n  n 

P™  (sin  ij.'  )  are  the  associated  Legendre  functions  defined  by 


Pm 

n 


(  X  }  a. 


2  m /Z 
(1  -  x  ) 


dx 


P  (x) 


In  practice,  the  program  is  capable  of  utilizing  all  zonal,  tes¬ 
seral,  and  sectorial  terms  thi  ough  (ZO,  ZO). 
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Drag  Considerations 


The  drag  effect  is  computed  from  a  double  interpolation  of  Jacchia's 
model  atmosphere  which  gives  density  as  a  function  of  height  and  exo¬ 
spheric  temperature.  The  table  contains  values  for  heights  ranging  between 
120  and  1000  kilometers  and  temperatures  between  2400°  K  and  600°  K, 

For  heights  above  1000  kilometers,  an  exponential  extrapolation  is  used. 

The  actual  density  calculation  proceeds  as  follows. 

The  exospheric  temperature  (T)  is  computed  from  Jacchia's  model 


for  the  diurnal  variation, 

T  =  T^  1 1  +  R  ^sin  6+  (cosmr|  -  sinm0)  cos*1  -j-  J  | 
where  T^  =  nightime  temperature 

ri  *  \  (  6  -  ) 

6  -  |  (6  +  6e) 

r  =  (a  -  a  J  +  c  +  p  sin  (i  -  (^  +  7  ) 

where  a,  6  *  right  ascension,  declination  of  the  satellite 
a  ,  6  -  right  ascension,  declination  of  the  sun 

Cl/  w 

R,  m,  n,  3.  p.  7  are  constants  determined  by  Jacchia  from  drag 
studies  using  the  Nicolet  II  model. 

The  height  of  the  satellite  is  computed  by  subtracting  the  radius  of 
the  earth  at  the  geocentric  latitude  of  the  subsatellite  point  from  the  geo¬ 
centric  distance  of  the  satellite.  The  radius  of  the  earth  is  approximated  by 


r  =  a  -  f  sin  6 
e 

vhere  6  *  geocentric  latitude 

f  *  21.  382  kilometers  (difference  between  equatorial  and  polar 

radii  of  earth) 

a  *  6378.  165  kilometers  (equatorial  radius  of  earth) 

Two  dimensional  quadratic  interpolation  in  the  table  is  then  used  10 


achieve  the  desired  density. 


Integration  Algorithm 

The  integration  algorithm  used  in  the  program  is  a  sixth  order 
predictor  -  corrector  scheme  which  is  initiated  by  a  classical  fourth 
order  Runge-Kutta  method.  The  algorithm  is  documented  separately 
in  Appendix  B  of  this  report. 

3.  Z  Differential  Orbit  Correction 


Each  observation  to  be  processed  gives  rise  to  £  or  3  (depending  upon 

3 

the  number  of  measured  quantities)  equations  of  condition  of  the  form  : 


cos  6  cos  a  cos  6  sin  a  sin  6 


\ 


3x1 


-  sin  a  cos  q 

P  P 

-  sin  6  cos  a  sin  3  sir>  a  cos  6 

3  x  3 


/  3x  3x 

3x  3y 
o  o 

3v  3y 

3x  3y 
o  o 


3z  3  z 

\3x  3y 
\  o  '  o 

3x7 


3A 

3A 

3z 

3A 


/M 

ldyo  1 

dz 

o 

dx 

o 

d*° 

dz 

1  ° 

7  x  1 


where 

Ap,  A  a,  A  6  are  the  differences  in  range,  right  ascension,  and 
declination  between  computed  and  observed  position. 

dx  ,  dy  ,  dz  ,  dx  ,  dy  ,  dz  ,  dA  are  the  desired  corrections  to 
o  o  o  o  o  o 

the  initial  position,  velocity,  and  atmospheric  model  in  order  to  minimize 
A  p  ,  Act,  and  A  6 . 

The  3x7  matrix  gives  the  partial  derivatives  of  present  position 
with  respect  to  initial  position,  velocity,  and  density.  These  derivatives 
are  computed  numerically. 

l'he  3x3  matrix  is  for  geometrical  purposes. 

When  all  observations  have  been  processed,  the  resultant  system 
of  equations  is  solved  in  the  least  squares  sense,  thereby  yielding  the 
desired  corrections.  Corrections  to  position  and  velocity  are  added 
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directly  to  initial  position  and  velocity,  while  the  correction  to  the  atmos¬ 
pheric  model  is  affected  by  modifying  each  interpolated  value  retrieved 
from  the  model.  The  entire  procedure  is  iterated  until  convergence  is 
achieved. 


Preliminary  Adjustments 


Prior  to  the  differential  correction  procedure,  a  preliminary 
adjustment  is  made  to  those  parameters  most  likely  to  contain  signifi¬ 
cant  errors,  namely,  mean  anomaly,  semi-major  axis,  and  the  atmos¬ 
pheric  model.  Given  an  initial  estimate  of  the  orbital  parameters  of  the 
vehicle,  the  numerical  integration  routine  is  called  to  generate  computed 
position  and  velocity  at  each  of  the  observation  times.  Differences  between 
computed  and  observed  position  and  velocity  are  transformed  into  mean 
anomaly  residuals^,  which  are  then  fitted  with  a  second  degree  polynomial. 
The  constant  term  is  used  as  a  correction  to  the  mean  anomaly,  the  coef¬ 
ficient  of  the  linear  term  is  used  as  a  correction  to  the  semi-major  axis, 
and  the  coefficient  of  the  quadratic  term  is  used  to  correct  the  atmospheric 
model.  (As  in  the  differential  correction  procedure,  the  correction  to  the 
model  is  affected  by  modifying  each  interpolated  value  retrieved  from  the 
model.  ) 
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Section  4 

Program  Construction 

The  program  contains  46  routines  each  of  which  is  labeled  and 
serialized  in  columns  73-80.  All  coding  is  in  FORTRAN  IV.  The 
name  and  a  brief  description  of  each  routine  is  shown  in  Table  4-1. 
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Table  4-  1 


Routines  in  the 

Name 

DRIVE 

ADMON 

RDPRCD 

STNEQX 

OBSET 

PREPOB 

SAO 

SPTRK1 

SPTRK2 

IPMAD 

POLYLS 

DOC 

PRTRES 

MARESD 

MATIV 

EQCDT 

GEOMTR 

PTLDRV 

INTCRL 

RUNKUT 

ABM 

RHSEQM 

MOON 


Cambridge  Atmospheric  Density  Numerical  Integration  Program. 

Function 

Activates  atmospheric  density  processing  routines. 

Main  control  routine. 

Reads  and  prints  data  cards. 

Reads  station  and  equinox  data. 

Groups  observations  into  sets. 

Performs  preliminary  calculations  on  observations. 

Reads  Smithsonian  format  observations. 

Reads  AFCRL  format  observations. 

Reads  SPADATS  format  observations. 

Performs  preliminary  adjustment  of  mean  anomaly, 
semi-major  axis,  and  atmospheric  model. 

Quadratic  curve  fitting  routine. 

Differential  correction  supervisor. 

Prints  residuals  after  differential  correction. 

Computes  mean  anomaly  residuals. 

Matrix  inversion  routine. 

Forms  equations  of  condition. 

Computes  geometry  matrix. 

Computes  partial  derivatives. 

Control  routine  for  numerical  integration. 

Runge-Kutta  integration  routine. 

Predictor -cor rector  integration  routine. 

Computes  right  hand  sides  of  equations  of  motion. 

Computes  lunar  position  (inactive  at  present). 
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Table  4-1  (continued) 


Name 

GEOPT 

LGNDR 

MODEL 

DENS 

DRAG 

SOL 

OSCTMN 

PVTELM 

MNTOSC 

ELMTPV 

SHPRDC 

XYZADR 

ADTAE 

TDIFR 

CLNTJD 

TINCR 

YRMNDY 

ND1NYR 

NDTOD 

TDCONV 

KEPLER 

ANGLEl 

ANGLE2 


Function 

Evaluates  geopotential  contribution. 

Computes  Legendre  polynomials  ana  associated  functions. 
Chooses  proper  atmospheric  model. 

Computes  density  from  Jacchia's  model  atmosphere. 
Computes  density  from  drag  data  (inactive  at  present). 
Computes  solar  position. 

Converts  position  and  velocity  to  mean  elements. 
Converts  position  and  velocity  to  osculating  elements. 
Converts  mean  elements  to  position  and  velocity. 
Converts  osculating  elements  to  position  and  velocity. 
Computes  short  periodic  perturbations. 

Converts  x,  y,  z  to  Gt,  6,  0 

Converts  a,  6  to  azimuth  and  elevation 

Computes  difference  between  two  times. 

Converts  calendar  date  to  modified  Julian  Date. 
Increments  a  given  time. 

Computes  date  from  day  of  year. 

Computes  number  of  days  in  a  year. 

Computes  day  of  year  from  date. 

Converts  a  date  to  output  format. 

Solves  Kepler's  equation. 

Reduces  an  angle  to  the  interval  (  -tr,+7r  ). 

Reduces  an  angle  to  trie  ii.t  rval  (  0,  Zn  ). 
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Section  5 

Summary  of  Investigations  and  Accomplishments 

This  section  summarizes  the  investigations  performed  and  the  accomp¬ 
lishments  achieved  during  this  contract. 

5.  1  Revised  Numerical  Integration  Algorithm 

A  revised  numerical  integration  algorithm  consisting  of  a  sixth  order 
predictor-corrector  method  which  is  started  by  a  classical  fourth  order 
Runge-Kutta  scheme  was  incorporated  into  the  program  resulting  in  a 
considerable  saving  of  computer  time.  See  Appendix  B  for  a  complete 
description  of  this  method. 

5.  2.  Elimination  of  Numerical  Integration  Re-Starts 

The  incorporation  of  the  predictor •  corrector  method  did  not  attain  the 
expected  reduction  in  computer  time.  It  was  soon  determined  that  this  was 
due  to  observational  data  in  the  time  interval  being  considered.  Each  ob¬ 
servation  being  processed  required  a  re-start  of  the  predictor-corrector 
method,  thereby  considerably  reducing  its  effectiveness.  Consequently, 
the  logic  of  the  numerical  integration  procedure  was  revised  so  as  to 
eliminate  the  Runge-Kutta  re-starts  at  each  observation.  This  modification 
reduced  the  running  time  to  its  expected  level. 

5.  3  Elimination  of  High  Order  Terms  from  Preliminary  Calculations 

Additional  logic  was  incorporated  into  the  program  which  allows  the 
user  to  eliminate  high  order  zonal,  tesseral,  and  sectorial  harmonics  for 
certain  preliminary  calculations  not  requiring  extreme  accuracy.  Specif¬ 
ically,  a  lower  order  geopotential  model  is  now  used  for  the  preliminary 
adjustment  procedure  (see  Section  3.3)  and  for  early  iterations  in  the 
differential  correction  procedure,  resulting  in  a  further  reduction  in 
running  time. 
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5.4  Analytic  Representation  of  the  Geopotential 

An  attempt  was  made  to  replace  the  existing  formulation  for  esti¬ 
mating  the  effects  of  the  geopotential  on  the  motion  of  a  satellite  by  an 
analytic  scheme.  The  successful  completion  of  this  task  would  result 
in  a  significant  reduction  in  the  running  time  of  the  program.  The  re¬ 
sults  of  this  investigation  are  presented  in  Section  6, 

5.  5  Calculation  of  Perigee  Data 

The  accuracy  of  the  program  was  improved  by  replacing  the  calcu¬ 
lation  of  the  final  perigee  time  and  height  with  a  more  appropriate  scheme. 
Originally,  the  final  perigee  time  and  height  were  obtained  from  the  final 
corrected  orbital  elements  by  setting  the  mean  anomaly  equal  to  zero  and 
by  using  the  keplerian  formula  for  the  period  for  computing  the  time  of 
perigee.  The  procedure  used  now  actually  integrates  from  the  time  of  the 
differential  correction  to  perigee  time.  Consequently,  the  effects  of  air 
drag  and  the  geopotential  are  now  taken  into  account.  The  perigee  time  is 
computed  by  an  iterative  approach. 

5.6  Recoding  of  Selected  Portions  of  the  Program 

Selected  portions  of  the  program,  including  the  entire  geopotential 
as  well  as  various  utility  routines  used  in  the  conversion  from  (to)  orbital 
elements  to  (from)  position  and  velocity,  were  recoded  in  double  precision. 
This  change  yielded  a  more  accurate  and  reliable  procedure. 

5.7  Inclusion  of  Drag  Model 

As  an  alternative  to  using  the  Jacchia  density  model,  a  drag  model 
was  incorporated  into  the  program,  allowing  the  user  to  supply  accelerometer 
type  data  as  a  function  of  time.  Thus  far  no  input  data  of  this  form  has  been 
made  available,  so  that  the  formulation  has  not  been  verified. 

5.  8  Modification  to  Output  Format 

The  output  format  of  the  program  was  modifed  significantly,  yielding 
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more  'readable'  results.  In  addition,  a  complete  listing  of  the  parameter 
table  (as  modified  by  the  user)  is  included  in  the  output.  This  listing  elim¬ 
inates  any  ambiguities  with  regard  to  initial  conditions. 

5,  9  Overlaying  of  Program 

A  significant  amount  of  core  storage  was  made  available  by  overlaying 
the  CADN1P  program  as  shown  in  figure  5-1. 

Link  0  is  the  main  link  and  resides  in  core  storage  at  all  times.  Links  1 
through  4  normally  reside  on  a  system  utility  device  and  are  read  into  a  re¬ 
served  area  of  core  storage  when  their  execution  is  required.  Similarly  links 
5  and  6  share  a  common  area  of  core  storage.  As  a  result  approximately 
4700  decimal  locations  have  been  saved.  These  locations  are  used  for  additional 
observation  buffer  space,  allowing  the  user  to  process  up  to  125  observations 
in  one  run. 
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Figure  5-1 


Structure  of  Overlayed  Version  of  CADN1P 


Link  (0) 
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Section  6 


Analytic  Representation  for  Geopotential 

In  the  present  version  of  the  CADNIP  Program,  the  variable  orbital 
parameters  are  linked  to  the  observations  by  means  of  numerical  integra¬ 
tion  of  the  rectilinear  equations  of  motion.  Such  a  procedure  has  some 
very  important  advantages.  Perturbations  due  to  air  drag,  which  are  quite 
difficult  to  treat  analytically,  are  handled  easily  with  numerical  integration. 
Numerical  integration  also  takes  interaction  and  higher-order  effects  into 
account  automatically.  It  is  extremely  time  consuming,  however,  and  the 
use  of  CADNIP  in  production  work  results  in  the  expenditure  of  large 
amounts  of  computer  time. 

Several  refinements  have  been  introduced  into  the  CADNIP  program 
(see  Section  5)  in  order  to  reduce  its  running  time.  The  basic  premise  of 
numerical  integration  has  never  been  altered,  however.  Simultaneous  and 
detailed  treatment  of  the  major  atmospheric  and  gravitational  effects  is  an 
essential  requirement  of  the  program.  Yet,  the  possibility  exists  that  the 
smaller  terms  in  the  gravitational  potential  could  be  treated  by  analytic 
methods.  Since  CADNIP  typically  utilizes  a  fairly  complete  geopotential 
model,  a  considerable  amount  of  computer  time  is  required  to  evaluate  the 
right  hand  sides  of  the  equations  of  motion  which  must  be  done  at  each  step 
of  the  integration.  In  analytic  form,  it  would,  of  course  be  necessary  to 
evaluate  the  perturbations  due  to  these  many  small  terms  only  at  the  times 
of  the  observations. 

Considerable  time  and  effort  has  gone  into  the  study  of  the  possibility 
of  an  analytical  representation  of  the  higher  order  terms  of  the  geopotential 
in  CADNIP.  This  work  has  centered  around  a  program  developed  to  compare 
the  perturbations  resulting  from  the  numerical  integration  scheme  used  in 
CADNIP  with  those  resulting  from  a  particular  analytic  development.  This 
program  determines  the  perturbations  under  numerical  integration  by  in¬ 
tegrating  over  a  specified  time  interval  with  a  specified  set  of  (low-order) 
harmonic  coefficients,  integrating  a  second  time  using  the  original  set  plus 
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any  specified  set  of  higher-order  coefficients,  and  then  taking  the  difference 
between  the  two  results  at  regular  intervals.  Differences  in  the  osculating 
Keplerian  elements  as  well  as  differences  in  position  are  evaluated  and  com¬ 
pared  with  those  computed  by  routines  based  on  the  analytic  development  of 

v  .  7,8 
Kaula. 


Kaula's  theory  has  been  successfully  applied  by  several  groups  working 
with  earth  satellite  observations.  It  is,  for  example,  used  in  the  DOI  pro¬ 
gram  of  the  Smithsonian  Astrophysical  Observatory  to  compute  the  effects 
of  the  tesseral  (and  sectorial)  harmonics  during  orbit  determinations  from 
precise  photographic  observations.  The  theory  will  not  be  described  in  com¬ 
plete  detail  here.  A  brief  sketch  of  the  development  will  be  given,  however. 

The  earth's  potential  field  is  usually  described  by  a  series  of  spherical 
harmonics  in  the  form: 
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where  r  is  the  radial  distance  from  the  center  of  mass,  <f>  the  latitude  and 

> .the  (east)  longitude;  the  J  >\  ,  or  C  ,  S  are  numerical  coefficients; 

nm  nm  nm  nm 

and  the  P  (simt)  are  the  Legendre  associated  functions  of  the  first  kind  de- 
nm 

fined  by: 
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Kaula  transforms  the  potential  from  spherical  coordinates  to  Keplerian  el¬ 
ements  ( <*» ,  (l,  i,  e,  M,  a)  with  the  result  that,  for  a  particular  term  of  the  po¬ 
tential: 
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nm 
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where  9  is  the  Greenwich  Sidereal  Time, 


S 

nmpq 


n-m  even 

cos  Q(n  -  2p)o>  +  (n  -  2p  +  q)M  +  m(ft  - 
n-m  odd 


n-m  even 

sin  £(n  -  2p)u  +  (n 
n-m  odd 


2p  +  q)M  +  m 


and  the  F(  i)  and  G(e)  are  rather  complicated  functions  which  are,  however, 
quite  familiar  in  classical  celestial  mechanics. 


With  the  potential  in  thie  form,  the  Lagrangian  equations  of  motion  can 

be  integrated  analytically  quite  easily  if  it  is  assumed  that  the  only  variations 

on  the  right  hand  sides  are  linear  variations  in  the  argument  of  S  .  The 

nmpq 

result  for  one  element  should  serve  as  ample  illustration.  For  a  particular 
sub-term,  vnmpq>  in  the  potential,  the  resulting  perturbation  in  the  ascending 
node,  say,  is  given  by: 


LU 


nmpq 


e  £  n+3,,  2. 

Ma  (1-e  )  sini 


(3F  /  3i)G  S 

nmp  npq  nmpq 


Q(n-2p)u  +  (n-2p+q)M  +  m(o  - 


where  S  is  the  integral  of  s  with  respect  to  its  argument.  The 

nmpq  nmpq 

number  of  AQ  necessary  to  form  An  is  limited  by  the  fact  that,  among 

nmpq  i  i  nm 

other  things,  G  %e  q  so  that  only  sub-terms  with  q  near  zero  need  be  con- 
npq 

sidered.  In  practice  a  subroutine  that  computes  and  tests  the  amplitudes  of 
a  sufficiently  large  range  of  sub-terms  and  flags  and  stores  those  that  are 
significant  within  a  stated  accuracy  requirement  is  used.  The  amplitudes  are 
normally  computed  only  once  (assuming  that  the  orbital  elements  are  known 
with  sufficient  accuracy)  for  the  time  interval  being  considered.  The  pertur¬ 
bations  at  any  time  are  then  obtained  by  multiplying  the  amplitudes  by  the  cor- 
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responding  S  or  S  computed  for  that  time.  This  procedure  is  ob- 

nmpq  nmpq  r 

viously  quite  economical  in  terms  of  computing  time. 

Certain  interactions  between  the  different  perturbations  have  not  been 
mentioned  that  slightly  complicate  matters.  Short  period  perturbations  in 
the  semi-major  axis  will  give  rise  to  additional  short  period  perturbations 
in  the  mean  anomaly  that  may  have  to  be  taken  into  accourt.  In  addition, 
the  perturbations  in  e  and  i  that  result  from  zonal  harmoiic  terms  interact 
with  the  secular  perturbations  in  ^  £1,  and  M  due  to  resulting  in  addi¬ 

tional  long  period  perturbations  in  u.n  >  and  M  that  may  be  important.  These 
effects  can,  fortunately,  be  included  very  easily  and  this  is  done  in  the  rou¬ 
tines. 

At  the  time  of  this  writing,  a  completely  satisfactory  comparison  be¬ 
tween  numerical  integrations  and  analytic  theory  has  not  been  obtained  and 
no  attempt  at  implementation  of  the  analytic  representation  in  CADNIP  has 
been  made.  While  the  two  methods  show  generally  good  agreement,  small 
differences  exist  which  thus  far  remain  unexplained.  Typical  results  ob¬ 
tained  from  the  two  methods  for  a  relatively  low  orbit  with  a  small  eccen¬ 
tricity  are  illustrated  in  figures  6-1  through  6  -  10.  The  first  6  figures 
show  the  effect  of  ^  on  the  Keplerian  elements  as  determined  by  the 
analytic  and  numeric  approaches.  The  next  3  figures  represent  the  position- 
a'  displacement  implied  by  the  two  methods  in  a  rectangular  coordinate 
system  whose  origin  is  at  the  center  of  the  vehicle.  Figure  6  -  10  is  a  mag¬ 
nification  of  the  differences  between  the  two  approaches  shown  in  figures 
6-7  through  6-9.  The  case  of  very  small  eccentricity  is,  of  course, 
quite  important  in  applications  of  CADNIP.  The  differences  decrease  for 
larger  eccentricity  but  do  not  entirely  disappear  until  eccentricities  on  the 
order  of  0.  1  are  reached.  These  computations  were  made  with  included 
in  the  reference  integration  used  in  determining  the  perturbations  under  nu¬ 
merical  integration.  An  interesting  result  is  that,  when  is  not  included, 

excellent  agreement  between  the  two  methods  seems  to  occur.  This  is  illus¬ 
trated  in  figures  6-11  through  6  -  20.  It  v<ould  appear,  then,  that  the  observed 
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differences  might  be  due  to  interactions  between  the  short-periodic  per¬ 
turbations  due  to  C^q  and  the  higher -order  terms  that  are  not  included  in 
the  analytic  theory.  The  short-periodic  perturbations  in  »,M  and  (rel¬ 
atively)  e  due  to  become  quite  large  for  small  eccentricity.  This 
hypothesis  was  tested  in  considerable  detail.  Such  interaction  terms  can 


be  computed  from: 


3s 

\0s.  =  y  "t  — ■  A1  s . 
2  1  (  1  1 


where  s.  is  any  one  of  the  Keplerian  elements  being  considered  with  re¬ 
spect  to  a  particular  sub-term  and  AjS,  t*le  interacting  perturbation  in 
any  one  of  the  elements.  The  3s./3s.'s  can  be  obtained  by  differentiation 
of  the  Lagrangian  equations  of  motion  after  substitution  of  the  various  de¬ 
rivatives  of  the  potential  with  respect  to  the  elements,  and  the  a  ^s.  can 
be  obtained  via  the  first-order  theory  outlined  above.  This  is  quite  simple 
in  theory  but  extremely  laborious  in  practice,  both  from  an  algebraic  and 
a  logical  point  of  view.  Nonetheless,  the  routines  were  modified  to  include 
the  interaction  terms,  in  both  directions,  between  the  short-periodic  per¬ 
turbations  due  to  C^0  and  any  particular  higher-order  term. 

The  result  of  this  work  was  that  the  interaction  terms  showed  rather 
small  amplitudes  (although  they  were  beginning  to  be  significant  in  terms 
of  the  most  accurate  observations)evan  for  low  eccentricity ,  andwerenotas 
large  as  the  observed  differences.  Moreover,  they  were  strictly  periodic 
ns  they  clearly  should  be.  The  observed  differences,  on  the  other  hand, 
are  somewhat  irregular  and  tend  to  grow  -with  time.  This  gives  some  rea¬ 
son  to  suspect  the  numerical  integration.  Careful  checking  and  comparison 
with  an  essentially  independent  formulation  have,  however,  failed  to  reveal 
any  fault  in  either  the  theory  or  coding  of  the  routines  involved  in  the  in¬ 
tegration. 

It  is  felt  that  the  possibility  of  analytic  representation  of  the  high-order 
in  the  geopotential  should  be  explored  further.  It  may  well  be  that 
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terms 


what  we  have  observed  is  the  beginning  of  failure  of  the  theory  in  the  case 
of  near  zero  eccentricity.  Such  failure  must  eventually  occur  when  the 
Keplerian  elements  are  used.  If  this  is  the  case,  the  problem,  which  would 
be  of  a  purely  mathematical  nature,  could  be  solved  by  a  proper  transfor¬ 
mation  of  the  orbital  elements  and  development  in  terms  of  the  Delaunay 
elements.  This  would,  clearly,  involve  a  considerable  amount  of  analysis 
and  programming.  Of  course,  if  this  were  determined  to  be  the  cause  of 
the  difficulty,  limits  might  be  established  and  the  existing  routines  incor¬ 
porated  in  CADNIP  for  use  with  larger  eccentricities.  Perhaps  consider— 
ation  should  be  given  to  both  of  these  approaches  in  future  work. 


-  20  - 


DEGREES*) .OE+5) 


Figure  6  -  1 


M  OMEGA  C  OMEGA  I  E  A 

0.000  0.000  0.000  H5.0000  .00500  7136.7028 

PERIOD  (MINUTES)  *  100.00 

RUN  FOR  C (2.2) 

(*)  NUMERIC 
(  +  )  ANALYTIC 
(O)  DIFFERENCE 


ij||| 


Figure  6-2 


OMEGA  C  OMEGA  I  E  A 

0.000  0.000  HS.0000  .00500  7136.7028 

PERIOD  (MINUTES)  =  100.00 

RUN  FOR  C(2.2) 

(*)  NUMERIC 
(  +  )  ANALYTIT 
(0.  DIFFrRENCE 


Figure  (>  -  3 


OMEGA  C  OMEGA  I  E  A 

0.000  0.000  45.0000  .00500  7136.7028 

PERIOD  (MINUTES)  ■  100.00 

RUN  FOR  C(2. 2) 

(1)  NUMERIC 
(  +  )  ANALYTIC 
(©)  DIFFERENCE 


I?'  i  I'Nfv  i'  !-t 1 i»T,s 1 ii ' ‘in' r|i  (fi  n|i  ifl'  ipf 


H 

0.000 


f 


Figure  6-4 

OMEGA  C  OMEGA  I 

0.000  0.000  HS. 0000 

PER  1 00  (MINUTES)  -  100.00 

RUN  FOR  C ( 2 • 2 ) 


E 

.00500 


(X)  NUMERIC 
t  +  )  ANALYTIC 
(O)  0IFrERENCE 


Figure 


OMEGA  C  OMEGA 

0.000  0.000 

PERIOD  (MINUTES 

RUN  FOR  C (2. 2) 

( #  >  NUMERIC 
l + )  ANALYTIC 
(©)  DIFFERENCE 


I 

MS. 0000 


TIME  (MINUTES) 


Figure  6-6 
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Figure  6-7 
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Figure  6-8 
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Figure  6-12 


OMEGA  C  OMEGA  1  E  A 

0.000  0.000  45.0000  .00500  7136.7028 

PERIOD  (MINUTES)  -  100.00 

RUN  FOR  C (2. 2) 

(I)  NUMERIC 
(  +  )  ANALYTIC 
(©)  DIFFERENCE 


J ( 2 )  ELIMINATED 


Figure  6-13 
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Figure  6-15 
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Section  7 


Recommendations 

The  program  developed  under  this  contract  has  proven  to  be  a 
valuable  tool  in  atmospheric  density  research.  The  improvements  and  re¬ 
finements  incorporated  in  CADNIP  have  resulted  in  a  significant  reduction 
in  computer  running  time  while  yielding  more  accurate  and  reliable  results. 
One  area  of  the  research  effort  has  been  disappointing  however.  The 
inability  to  successfully  incorporate  in  CADNIP  the  analytic  formulation 
for  the  geopotentiai  has  been  a  source  of  anguish  and  frustration.  Con¬ 
sidering  the  time  and  effort  already  devoted  to  this,  and  considering  the 
savings  in  computer  time  that  would  be  realized  if  this  were  successful, 
it  would  seem  that  additional  effort  in  this  direction  as  outlined  in  section  6 
is  warranted. 

An  additional  area  of  development  which  is  recommended  is  the 
inclusion  of  more  recent  and  sophisticated  atmospheric  models  in  the  pro¬ 
gram.  This  would  be  relatively  easy  to  accomplish  since  the  mechanism 
or  logic  for  the  insertion  of  new  models  was  essentially  incorporated  in 
CADNIP  when  the  drag  model  described  in  section  5.  7  was  introduced. 
Alternatively,  the  CADNIP  program  may  be  used  as  a  test  vehicle  for  at¬ 
mospheric  models.  Given  a  set  of  observations,  the  program  would  be 
run  using  a  number  of  different  models.  The  model  which  yields  an 
epbecneris  closest  to  known  positions  would  be  considered  most  repre¬ 
sentative  of  the  atmosphere. 

In  the  present  version  of  the  CADNIP  program,  the  processing 
of  a  set  of  observations  is  terminated  by  the  completion  of  the  differential 
correction  procedure.  It  is  possible  however  under  certain  circumstances 
to  deduce  additional  density  information  from  the  observations.  A  plot  of 
the  mean  anomy  residuals  (which  are  expected  to  appear  as  random  noise) 
may  reveal  trends  which  could  yield  additional  density  data.  Thus  it  is 
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suggested  that  a  post-processor  be  included  in  CADNIP  which  will 
examine  the  final  residuals  in  an  attempt  to  isolate  any  previously  ignor 
data  such  as  time  dependent  density  information. 
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This  section  of  the  report  is  concerned  with  the  actual  use  of  the 
program  and  includes  descriptions  of  all  input  data  as  well  as  operating 
instructions . 

The  program  requires  four  types  of  data  for  every  run:  satellite 
observations;  station  coordinates  and  related  data;  a  model  atmosphere; 
and  control  cards  which  specify  the  parameters  of  a  run. 

A.l  Observation  Preparation 

The  present  version  of  the  program  is  capable  of  processing  ob¬ 
servational  data  prepared  in  one  of  three  formats;  Smithsonian  format; 
AFCRL  format;  or  SPADATS  format. 

A.  1.  1  Smithsonian  Format 

Each  observation  is  punched  onto  a  card  which  contains  the  follow¬ 
ing  information.  (This  description  will  be  limited  to  the  portion  of  the 
card  interpreted  by  the  program.  A  complete  description  may  be  found 
in  reference  9  .  ) 


Field  Columns  Contents  -  Description 


1 

14-17 

Station  number. 

2 

18-23 

Year,  month  and  day  of  the  observation 
with  two  columns  for  each. 

3 

c4-28 

Hour  of  observation. 

4 

26-27 

Minute  of  observation. 

5 

28-33 

Seconds  of  observation.  The  decimal  point 
is  assumed  to  be  between  columns  29  and  30 

V 

24-  36 

Degrees  portion  of  azimuth  or  hours  por¬ 
tion  of  right  ascension. 

7 

37-38 

Minutes  portion  of  azimuth  or  right  ascen- 
s  ion. 

8 

39-43 

Seconds  portion  of  umulh  or  right  ascen¬ 
sion,  The  decimal  point  is  assumed  to  be 
between  columns  40  and  41. 
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Field 

Columns 

Contents  -  Description 

9 

44-46 

Degrees  pottion  of  elevation  or  declination. 

10 

47-48 

Minutes  portion  of  elevation  or  declination. 

1  1 

49-52 

Seconds  portion  of  elevation  or  declination. 
The  decimal  point  is  assumed  to  be  be¬ 
tween  columns  50  and  51. 

12 

56 

Type  of  observation 

0^  right  ascension  and  declination 

1  or  3,  azimuth,  elevation  and  possibly 
range 

1  3 

57 

Code  for  equinox  of  reference 

1,  equinox  of  1855.  0 

2,  equinox  of  1875.0 

3,  equinox  of  1900.0 

4,  equinox  of  1950.0 

14 

59-63 

Range  in  kilometers.  The  decimal  point 

is  assumed  to  be  between  columns  59  and 
60. 


15  64  Power  of  ten  by  ,  'hich  to  multiply  range 

field  to  obtain  true  range. 

The  following  comments  will  be  useful  to  the  reader  in  preparing  a 
Smithsonian  format  observation  tape. 

Each  observation  is  read  according  to  the  following  FORMAT  statement. 


FORMAT  ( 13 X , 1 4 , F6 . 0 , 2F2 . 0 , F6 . 4 ,F3 . 0 , F2 . 0 , F5 . 3 , F3 . 0 , F2 . 0 , F4 . 2 , 3X , 21 1 ,  IX , F5 . 4 , 1 1 ) 


When  all  desired  observations  have  been  prepared  on  cards,  the 
cards  are  pul  on  a  tape  in  ascending  chronological  sequence  with  a  blank 
c  ard  following  the  last  observation. 

A  prefix  of  1900  is  assumed  for  the  year  field  of  each  observation. 

The  equinox  reference  code  (field  13)  is  only  interpreted  for  right 
asiersion,  declination  type  observations,  while  the-  range  field  is  only  inter¬ 
preted  for  azimuth,  elevation  type  observations.  An  azimuth,  elevation 
type  observation  is  converted  to  a  range,  azimuth,  elevation  type  observa¬ 


tion  it  the  range  field  is  greater  than  zero. 


Since  SAO  uses  the  same  format  for  field  reduced  and  photo- reduced 


observations,  it  is  possible  to  process  either  type  observation,  but  not 
both  together.  This  is  true  since  different  time  systems  are  used  for 
field  reduced  and  photo-reduced  observations.  Since  the  program  does 
not  determine  the  type  of  instrumentation  used  in  the  observation,  the  user 
is  required  to  separate  the  types  manually.  This  can  be  done  by  sorting 
the  observation  cards  on  columns  8  which  will  contain  a  7  for  all  photo- 
reduced  observations. 

A.  1.2  AFCRL  Format 

Each  observation  is  punched  onto  a  card  which  contains  the  follow¬ 
ing  information.  (This  description  will  also  be  limited  to  the  portion  of 
the  card  interpreted  by  the  program.  ) 

Field  Columns  Contents-Description 


1 

9-12 

Station  number. 

2 

14 

Units  digit  of  year  of  observation 

3 

16-17 

Month  of  observation. 

4 

19-20 

Day  of  observation. 

5 

22-23 

Hour  of  observation. 

6 

25-26 

Minute  of  observation. 

7 

28-32 

Seconds  of  observation.  The  decimal 
point  is  assumed  to  be  between  columns 

29  and  30. 

>  8 

34-39 

Elevation  in  degrees.  The  decimal 
point  is  assumed  to  be  between  columns 

36  and  37. 

9 

41  -46 

Azimuth  in  degrees.  The  decimal 
point  is  assumed  to  be  between  columns 

43  and  44. 

10 

48-56 

Range  in  kilometers.  The  decimal 
point  is  assumed  to  be  between  columns 

53  and  54. 

l _ 
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1 

1 

4 

The  following  comments  will  be  useful  to  the  reader  in  preparing 
an  AFCRL  format  observation  tape. 

Each  observation  is  read  according  to  the  following  FORMAT 
statement: 

FORMAT  (8X,  14,  IX.  FI.  0,  4(1X,  F2.0),  IX,  F5.  3,  2(1X,  F6.3),  IX,  F9.  3) 
When  all  desired  observations  have  been  prepared  on  cards,  the 
cards  are  put  on  a  tape  in  ascending  chronological  sequence  with  a  blank 
card  following  the  last  observation. 

A  prefix  of  I960  is  assumed  for  the  year  field  of  the  observation. 

A.  1.3  SPADATS  Format 

Each  observation  is  punched  onto  a  card  which  contains  the  follow¬ 
ing  information.  (This  description  will  also  be  limited  to  the  portion  of 
the  card  interpreted  by  the  program.  A  complete  description  may  be  found 
in  reference  10.  ) 

PM e  1  d  Columns  Contents  -Description 

1  1  Classification  indicator. 

S=Secret,  C=Confidential,  U=llnclassified. 

2  7-9  Station  number. 

3  10-1  1  Last  two  digits  of  year  of  observation. 

4  12-14  Day  of  year  of  observation. 

5  15-16  Hour  of  observation. 

b  17-18  Minute  of  observation. 

7  19-20  Second  of  observation. 

8  21-23  Thousandths  of  a  second. 

9  24  Tens  digit  of  ele vation/declination  field  in  degrees. 

This  field  may  also  be  overpunched  with  a  minrt'  punch. 

10  25-29  Remainder  of  elevation/ declination  field  in  degrees, 

11  31-32  Hundreds  and  tens  digit  of  azimuth  field  or  hours 

portion  of  right  ascension  field.  Azimuth  measure¬ 
ments  are  in  degrees, 

12  33-  34  Units  and  tenths  digits  of  azimuth  field  or  minutes 

portion  of  right  ascension  field. 

15  35-37  Hundredths,  thousandths,  and  ten  thousandths  digits 

of  azimuth  field  or  tens,  units  and  tenths  digits  o£ 
seconds  portion  of  right  ascension  field. 
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Field 


Columns 


Contents  -  Description 


Range  field  in  kilometers.  The  decimal  point  is 
assumed  to  be  between  columns  40  and  4L 

Power  of  ten  by  which  to  multiply  field  14  to 
obtain  true  range. 

Type  of  observation 
1  -  azimuth,  elevation 
2,3,4  -  range,  azimuth,  elevation 
5  -  right  ascension,  declination 
0,  6,  7,  8,  9  -  not  recognized 

Equinox  indicator  for  right  ascension,  declination 
observations 

1  -  Equinox  of  1900.0 

2  -  Equinox  of  1920.0 

3  -  Equinox  of  1950.0 

4  -  Equinox  of  1975.0 

5  -  Equinox  of  2000.  0 

6  -  Equinox  of  1850.0 

7  -  Equinox  of  1855.  0 

8  -  Equinox  of  1875.0 

9  -  Equinox  of  1960.  0 
0  -  Not  recognized 

The  following  comments  will  be  useful  to  the  reader  in  preparing  a 
SPADATS  format  observation  tape. 

Each  observation  is  read  according  to  the  following  FORMAT  statement. 

FORMAT  (A1,5X, 13 ,12 , 13 , 2F2 .0 ,F5 . 3 ,A1 ,F5 . 4 , IX , 2F2 . 0 ,F3 . 1 , IX, F7 . 5 , 1 1  , 2RX, 211) 

When  all  desired  observations  have  been  prepared  on  cards,  the  cards 
are  put  on  a  tape  in  ascending  chronological  sequence  with  a  blank  card 
following  the  last  observation. 

A  prefix  of  1900  ie  assumed  for  the  year  field  of  the  observation.  < 

A.  2  Station  and  Equinox  Data  Preparation  I 

The  station  data  tape  is  generated  by  an  off-line  card-to-tape  opera-  J 

tion.  The  tape  includes  three  sections  of  data,  each  cf  which  must  be 

I 

followed  by  a  blank  card.  | 

f 

Section  one  consists  of  a  list  of  classified  stations  along  with  their  j 

unclassified  codes  which  are  used  in  all  output  which  references  these 
stations.  Each  card  of  section  one  contains  two  numbers,  where  the  first 


1  4 

39-45 

15 

46 

16 

75 

17 

76 
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number  is  the  classified  number  of  a  station,  and  the  second  number 
is  its  unclassified  code.  This  type  card  is  punched  in  accordance  with 
the  following  FORMAT  statement: 

FORMAT  (215) 

Notice  that  if  no  classified  station  data  exists,  the  first  record  on  the 
station  tape  will  be  a  blank  card. 

The  second  section  of  data  consists  of  one  card  for  each  station 


tor  which  there  is  an  observation.  Each  card  contains  the  station  number, 
latitude  (in  degrees),  longitude  (in  degrees),  and  height  above  sea  level 
(in  meters)  which  are  punched  in  the  above  specified  order  in  accordance 
with  the  following  FORMAT  statement: 

FORMAT  (14,  F7.4,  F8.  4,  F6.  0) 

Latitudes  may  range  from  -90  to  +90  with  the  northern  hemisphere 

assuming  positive  latitude  values.  Longitudes  may  range  from  -360  to 

+  360  with  longitudes  which  are  west  of  Greenwich  being  considered  positive. 

The  format  for  station  data  cards  is  in  full  agreement  with  the 

station  deck  which  may  be  obtained  from  Space  Track.  The  distributed 

Smithsonian  station  data  deck  is  in  a  different  format,  however,  and  a 

program  is  required  to  make  the  necessary  conversion.  This  program 

has  been  written  and  is  available  along  with  a  "converted"  Smithsonian 

station  deck  which  was  obtained  from  this  program. 

The  third  section  of  the  station  tape  contains  equinox  correction 

data  which  are  used  in  updating  right  ascension  (a)  and  declination  ('•) 

measurements  from  a  given  equinox  to  a  common  equinox.  An  equinox 

4 

correction  is  computed  from  the  following  formulas. 

Let 

At 

A  =  a  +  — —  (C^  +  tan  6  sin  a  ) 

and 

D  =  -  +  cos  a 

Then  the  updated  values  for  i  and  -  are  giv^n  by 


where 


a  =  the  difference  in  time  between  the  desired  equinox  and  a 
given  equinox 

Cf,  C^>  C3.  C4  are  the  appropriate  precessional  constants. 

The  information  associated  with  each  given  equinox  is  punched 
on  a  card  which  contains  the  year  of  the  equinox,  C^,  C^,  C^.  and 
in  the  specified  order,  in  accordance  with  the  following  FORMAT  statement: 

FORMAT  (5F12.0) 

The  order  of  the  equinox  cards  is  important  since  the  routines 
which  read  observational  data  refer  to  the  equinox  cards  by  number  rather 
than  by  year  of  equinox.  Consequently,  the  present  version  of  the  program 
requires  that  the  first  nine  equinox  cards  on  the  station  tape  correspond 
to  the  equinoxes  of  1355,  1875,  1900,  1950,  1850,  1920,  I960,  1975,  2000. 
These  nine  cards  presently  contain  the  following  values: 


1 

Year 
1855. 0 

C1 

. 0127979 

C2 

.0055696 

S 

. 0128022 

C4 

. 0055683 

2 

1  875.  0 

.0127995 

.0055691 

. 0128029 

. 0055681 

3 

1900. 0 

.0128014 

. 0055686 

.0128039 

. 0055678 

4 

1950. 0 

. 0128053 

. 0055674 

.0128058 

. 0055672 

5 

1850. 0 

. 0127975 

. 0055697 

. 0128020 

. 0055684 

6 

1920. 0 

. 0128029 

. 0055681 

.0128047 

. 0055676 

7 

i960. 0 

.0128060 

. 0055671 

.0128062 

. 0055671 

8 

1975.  0 

.0128072 

. 0055668 

.0128068 

. 0055669 

9 

2000. 0 

.0128091 

. 0055662 

.0128078 

. 0055666 

Because  of  storage  limitations,  there  is  a  maximum  number  of 
cards  which  may  be  present  in  each  of  the  three  sections  of  the  station 
tape.  The  present  version  of  the  program  allows  for  twenty-five  classi¬ 
fied  stations  in  section  one,  one  hundred  stations  in  section  two,  and 
ten  equinoxes  in  section  three.  In  the  event  a  station  list  contains  more 
than  the  allowable  number  of  entries,  it  will  be  necessary  to  isolate  the 
portion  of  the  station  list  referred  to  by  the  observation  tape.  For  this 
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purpose,  a  preprocessing  program  has  been  written  which  reads  observa¬ 
tions  and  lists  all  stations  referred  to  by  these  observations.  (This  program 
performs  two  additional  functions.  The  observation  tape  is  checked  for 
chronological  correctness  and  for  the  validity  of  numeric  fields.  ) 

A.  3  Density  Model  Preparation 

The  program  currently  uses  Jacchia’s  model  atmosphere^  This 
model  atmosphere  is  contained  on  534  cards  which  give  densities  for 
heights  between  120  and  1000  kilometers  in  steps  of  10  kilometers  and  for 
temperatures  between  2400  and  600  degrees  Kelvin  in  steps  of  50  degrees. 
The  routine  which  reads  this  model  (DENS)  has  been  made  quite  flexible 
by  the  addition  of  one  card  at  the  beginning  of  the  model  which  lists  various 
parameters  of  the  model.  This  capability  allows  future  expansions  or  con¬ 
tractions  of  the  model  without  requiring  any  revisions  to  the  DENS 
routine.  The  additional  card  contains  the  following  information. 


Field 


1 

2 

3 

4 

5 

6 

7 

8 


Value 

37 

2400 

-50 

100 

89 

120 

10 

50 


De  scription 

Number  of  temperature  values. 

Initial  temperature  value. 

Increment  in  temperature  entries. 

Allowable  extrapolation  (below  end  of  model) 
for  low  temperatures. 

Number  of  height  values. 

Initial  height  value. 

Increment  in  height  entries. 

Factor  used  in  high  height  extrapolation. 


This  card  is  read  according  to  the  following  FORMAT  statement. 

FORMAT  (  2  (16,  3F6.  0)  ) 

The  model  itself  contains  seven  entries  per  card  with  density  values 
lor  all  temperatures  at  a  given  height  requiring  six  cards.  (The  sixth 
card  contains  only  two  entries.)  The  density  values  are  recorded  as  com¬ 
mon  logarithms. 
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The  model  is  read  according  to  the  following  FORMAT  statement. 
FORMAT  (7F1  0.  0) 

The  density  tape  required  by  the  program  is  prepared  by  placing 
on  a  tape  the  special  card  (described  above)  followed  by  the  534  cards  of 
the  model. 

A.  4  Specification  Cards 

This  section  describes  the  method  of  preparing  the  cards  which 
specify  the  parameters  of  a  run. 

There  are  four  types  of  cards  which  may  be  used.  These  are  now 
described. 

A.  4.  1  Remark  Card 

This  type  of  card  is  used  to  introduce  any  desired  information  into 
the  output  of  a  run.  Remark  cards  if  present  are  the  first  cards  of  a  run. 

A  remark  card  is  identified  by  the  word  REMARK  in  columns  1-6,  with 
the  remainder  of  the  card  being  arbitrary. 

A  .  4.  2  Starting  Elements  Cards 

Following  any  remark  cards  is  a  set  of  two  cards  which  contains  the 
initial  or  estimated  orbital  parameters  of  the  satellite  at  a  specified  time. 
The  specified  time  should  be  near  the  beginning  of  the  period  of  time  to 
be  processed.  The  first  card  of  this  set  is  prepared  as  follows: 

Field  Columns  Contents  -Description 


1 

18-19 

Two  low  order  digits  of  year  of  elements. 

l 

21-22 

Month  of  elements. 

3 

24-25 

Day  of  elements. 

4 

27-28 

Hour  of  elements. 

* 

5 

30-31 

Minute  of  elements. 

:  6 

3  3-38 

Second  of  elements. 

7 

5  3-  58 

Nighttime  exospheric  temperature  at  time 
specified  in  fields  1-6. 

8 

65-72 

2 

Area/mass  ratio  (cm  /  g) 
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1 

1 

_ s 

13-24 
25-36 
37-48 
49-60 
61  -72 


Mean  anomaly  in  degrees. 

Argument  of  perigee  in  degrees. 

Right  ascension  of  ascending  code  in  degrees. 
Inclination  angle  in  degrees. 

Eccentricity. 

Semi-major  axis  in  kilometers. 


All  angular  elements  should  be  reduced  to  the  interval  (0,  360)  with 
the  additional  restriction  that  the  inclination  angle  be  limited  to  the  inter¬ 
val  (0,  180). 

This  card  is  read  according  to  the  following  FORMAT  statement. 

FORMAT  (4F12.  4.  F12.7.  F12.  3) 

A  typical  set  of  starting  element  cards  might  be 
TIME  OF  ELEMENTS  65  01  02  0  3  04  05.  678  NIGHT  TEMP.  678.  A/M.  0123 


123. 45 


123.  45 


123.  45 


12. 345 


.00345  6543.2 


Notice  that  the  first  card  is  made  quite  ''readable1  by  the  addition 
of  text  in  uninterpreted  portions  of  he  card. 

A. 4,  3  Change  Card 

At  the  beginning  of  each  run,  the  program  assigns  values  to  a 
number  of  parameters  which  are  not  likely  to  change  from  run  to  run.  The 
value  of  any  of  these  parameters  may  be  modified  by  the  user  by  means  of  a 
change  card.  The  contents  of  the  card  are  as  follows: 


!•'  i  c  1  d 


Columns 


Content  s  -  Di  s c  r  lptiun 


CHANGE 


2  4-37 


Number  ol  parameter  to  be  modified. 
Indication  of  parameter  mode. 
Modified  value  of  the  parameter. 


A 


If  field  three  is  blank,  the  program  assumes  that  the  parameter 
to  be  changed  is  a  real  variable;  otherwise  the  program  assumes  that  the 
parameter  to  be  changed  is  an  integer  variable. 

The  change  card  is  read  according  to  the  following  FORMAT  statement. 

FORMAT  <A6,  10X,  13,  Al.  3X,  EJ4.  7) 

Notice  that  field  four  must  be  in  the  real  mode  regardless  of  the  mode 
of  the  parameter  to  be  changed. 

A  typical  change  card  might  be; 

CHANGE  PARAMETER  3*  TO  10.  0E  +  00 

which  informs  the  program  that  the  observation  tape  to  be  processed  is  to 
be  found  on  logical  tape  number  10. 

A  complete  description  of  the  parameter  array  can  be  found  in 
paragraph  A.  5  of  this  appendix. 

Change  cards  if  present  come  immediately  after  the  starting  element  card 
A  ,  4.  4  Run  Ca  rd 

Following  any  change  cards  is  a  run  card  which  is  the  last  card  of  each 
run.  The  run  card  is  prepared  as  follows; 


Field 

1 

2 

3 

4 

5 

6 

7 

8 
9 

1  0 


Columns 

10-11 

13-14 

16-17  Start  time  of  processing 
19-20 
22-23 
28-29 
31-32 

34-35  End  time  of  processing 

37-38 

40-41 


Contents -Description 

Two  low  order  digits  of  year 

Month 

Day 


Hour 
\,Minute 

fTwo  low  order  digits  of  year 
Month 
Day 
|  Hour 
, Minute 


: 
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Field 


Columns 


Contents  -Desc  ription 


1  1  56  Observation  format  indicator 

1,  Smithsonian  format 

2,  AFCRL  format 

3,  SPADATS  format 

13  73  Atmospheric  model  indicator 

0,  1,  3-9,  Jacchia  model 
3,  Drag  model  (inactive) 

The  run  card  is  read  according  to  the  following  FORMAT  statement. 

FORMAT  (8X ,  5F3.C,  3X,  5F3.0,  13X,  I  3,  14X,  13) 

A  typical  run  card  might  be; 

RUN  FROM  69  01  02  00  00  TO  69  01  03  00  00  FOR  OB.  TYPE  3  AND  MOD.  TYPF.  1 

which  re  ^csts  the  program  to  process  all  observations  on  the  tape  con'aining 
SPADATS  format  observations  which  were  recorded  on  January  3,  1969 
using  the  Jacchia  model  atmosphere. 

A  .4.6  End  C,i  rd 

The  la.' :  run  in  a  sequence  of  CADNIP  runs  must  be  followed  by  a  card 
containing  the  letters  END  in  rolumns  J-3  and  blanks  in  columns  4-6. 

A.  6  Description  of  Parameter  Table 

The  parameter  table  contains  parameters  whose  values  are  subject  to 
change  by  the  user  by  means  of  the  change  card  as  described  in  paragraph  A.  4.  3 
of  'his  appendix.  As  the  parameter  table  iscompletely  initialized  by  the  program 
at  the  beginning  of  each  run,  the  effect  of  CHANGE  cards  is  not  carried  over  from 
run  to  run. 

Since  the  user  is  required  to  know  the  mode  of  any  parameter  he  desires 
to  modify,  the  following  convention  is  adopted.  A  parameter  is  a  real  parameter 
it  and  only  if  its  value  is  printed  in  this  report  as  a  number  containing  a  decimal 
point . 
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Parameter 

Assigned 

Number  _ 

Value 

Parameter  Use 

1 

1 

Logical  number  of  station  aid  equi  iox  tape. 

2 

2 

Logical  number  of  atmospheric  model  tape. 

3 

3 

Logical  numberof  observation  tape. 

6 

6378 .  165 

Equatorial  radiu3  of  earth  in  kilometers. 

7 

1. 0/298. 3 

Flattening  of  earth. 

8 

806.  8136 

Canonical  unit  of  time  in  seconds. 

9-10 

690101 . 0, 0. 0 

The  date  1  Tan.  1969  in  internal  format. 

11 

. 279073160 

Sideral  angle  of  Greenwich  at  1  Jan.  1969 
in  revolutions. 

12 

.  2 7 '5  79  093E-2 

Modified  value  for  rotation  rate  of  Green-wich 
in  revolutions  per  day. 

14-20 

0 . 0 

Lunar  orbital  elements  at  a  specified  time. 

21 

0 . 0 

Ratio  of  mass  of  moon  to  mass  of  earth. 

22 

2 . 2 

Coefficent  of  drag  (C_). 

23 

1  .0 

Ratio  of  rotation  rate  of  atmosphere  to  ro¬ 
tation  rate  of  earth. 

24 

l  .  0 

Initial  factor  to  apply  to  atmospheric  model. 

29 

0  .  1 

Smallest  allowable  factor  to  apply  to  model. 

10.0 

Largest  allowable  factor  to  apply  to  model. 

27 

.  28 

r" 

28 

1  .  5 

m 

29 

2.  5 

n 

6  constants  appearing  in  Jacchla' 

30 

-45 . 0 

6 

*  model  for  the  diurnal  variation. 

31 

12 . 0 

P 

32 

33 

34 

45 . 0 

0 . 0 

0.  0 

y  - 

] 

Factors  controlling  assumed 
*  latitude  of  the  center  of  the 

diurnal  atmospheric  bilge, 

35 

0.125 

Desired  accuracy  of  9olar  position  (in  days). 

37 

15.0 

Runge-Kutta  time  step  in  seconds  for  regular 
integration 

38 

4 

Ratio  of  predi ctov -cor rector  time  step  to 
Runge-Kutta  time  step  for  regular  integration. 

49 

0.  375 

Desired  time  span  in  days  of  a  sei  of  obser¬ 
vations  . 

A  -1  3 


Pa  ramete  r 
Xuinbe  r 


Assigned 

Value 


Parameter  Use 


41 

0.  125 

Amount  of  da /s  by  which  to  modify  time 
span  of  a  set  of  observations  if  necess  iry, 

42 

0.  25 

Minimum  allowable  time  span  in  days  of  a 
set  of  observations. 

4  1 

2.0 

Maximum  allowable  time  span  in  days  of 
a  set  of  observations. 

46 

6 

Minimum  number  of  observations  necessiry 
in  preliminary  adjustment  procedure. 

47 

1 

Number  of  iterations  in  preliminary  adjust¬ 
ment  procedure. 

43 

3.  0 

The  product  of  this  parameter  and  sigma 
determine  tolerance  used  in  rejection  of 
observations  during  preliminary  adjust¬ 
ment  procedure. 

5  ' 

7 

Minimum  number  of  observations  necessary 
in  differential  correction  procedure. 

52 

0 

Parameter  controlling  intermediate  output 
from  differential,  correction  procedure. 

5  ) 

14 

Minimum  number  of  equations  of  condition 
in  differential  correction  procedure. 

54 

10 

Maximum  allowable  number  of  iterations 
in  differentia1  correction  procedure. 

55 

> 

Number  of  iterations  in  differential  cor¬ 
rection  procedure  in  which  to  compute 
pa  -  ti  a  1  de i  va :  i  ves . 

5  o 

15.  0 

R.K.  time  step  in  seconds  for  computing 
par'ia1  derivatives  of  present  position  with 
respect  to  initial  position  and  velocity. 

57 

l  5.  0 

R.K.  time  step  in  seconds  for  computing  par¬ 
tial  derivatives  of  p resent  position  with  re  ¬ 
spect  to  density. 

53 

4 

Ratio  of  predictor -corrector  time  step  to  R.K. 
time  step  for  computing  partial  derivatives  of 
present  position  with  respect  to  initial  position 
and  velocity. 

54 

4 

Ratio  of  predictor -cor rector  time  step  to  R.K. 

time  step  for  computing  par'ial  derivatives  of 
present  position  with  respect  to  density. 
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Parameter 

Number 

Assigned 

Value 

Parameter  Use 

61 

1. 0E-  6 

Increment  used  in  estimating  partial  deri¬ 
vatives  of  present  position  with  respect  to 
initial  position  and  velocity. 

62 

1.  OE-  3 

Increment  used  in  estimating  partial  deri¬ 
vatives  of  present  position  with  respect  to 
density. 

65 

10 

Number  of  iterations  in  differential  correction 
procedure  in  which  to  reject  bad  observations. 

66-67 

5.  0 

Tolerances  in  degrees  used  for  rejection  of 
right  aecension/azimuth  and  declination/ele¬ 
vation  measurements  during  first  iteration 
of  differential  correction  procedure. 

68 

500.0 

Tolerance  in  kilometers  for  rejection  of  range 
measurements  during  first  iteration  of  differ¬ 
ential  correction  procedure. 

69 

50.  0 

Tolerance  in  kilometers  for  rejection  of  range 
measurements  during  second  and  subsequent 
iterations  of  differential  correction  procedure. 

70 

3.0 

The  product  of  this  parameter  and  sigma  de¬ 
termines  tolerance  used  in  rejection  of  right 
ascension/azimuth  and  declination/elevation 
measurements  during  second  and  subsequent 
iterations  of  the  differential  correction  pro¬ 
cedure. 

71 

5.0 

Tolerance  used  to  determine  whether  differential 
correction  procedure  has  diverged. 

72 

.  01 

Tolerance  used  to  determine  desired  convergence 
in  differential  correction  procedure. 

72 

.  1 

Tolerance  used  to  determine  acceptable  con¬ 
vergence  in  differential  correction  procedure. 

74 

15.  0 

Tolerance  in  minutes  of  arc  used  to  deter¬ 
mine  acceptability  of  the  differential  cor¬ 
rection  results. 

75 

.  1 

Tolerance  used  to  determine  acceptability 
of  the  final  density. 

76 

xxxx. 

Equinox  to  which  observations  are  updated. 

77 

XXX. 

Standard  height  density. 
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Parameter 
Numbe  r 


Assigned 

Value 


Parameter  Use 


79 


4 


Order  of  gravity  model  used  in  preliminary 
calculations. 


80 


8 


Order  of  gravity  model  used  in  calculations 
requiring  full  accuracy. 


r~ 

81 

8Z-83 
84  ; 

I 

i 

83-87  |  See  Table  A  -  1 

i 

l 

88  ; 


89-308' 


300-51 
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•ial  Harmonics 


Parameters  9-12  allow  the  user  to  specify  an  initial  position  and  a  rota¬ 
tion  rate  for  Greenwich  at  a  given  time.  Parameter  9  contains  a  lime  in  the 
form  of  DAY  +  100  X  MONTH  +  10000  X  YEAR.  Parameter  10  contains  a 
fraction  of  a  day.  Parameter  11  specifies  the  sidereal  angle  of  Greenwich 
at  the  time  specified  by  parameters  9  and  10.  Parameter  12  is  a  modified 
value  for  the  rotation  rate  of  Greenwich.  The  modified  value  which  is  es- 
senlially  @  -  1.0,  allows  for  the  retention  of  sufficient  accuracy  without 
resorting  to  double  precision  arithmetic. 

Lunar  effects  are  currently  ignored  by  the  program.  Consequently, 
parameters  14-21  are  set  to  zero. 

Parameter  24  will  be  modified  by  the  preliminary  adjustment  procedure 
and  the  differential  correction  procedure.  Parameters  25  and  26  act  as  lower 
and  upper  bounds  on  parameter  24.  This  feature  is  useful  in  controlling  com¬ 
putational  difficulties  which  may  arise  because  of  grossly  inaccurate  starting 
conditions  or  insufficient  drag  effects. 

In  the  special  case  where  parameters  25  and  26  are  both  set  equal  to  1.0, 
the  preliminary  adjustment  procedure  will  correct  the  mean  anomaly  and 
semi-major  axis  only,  while  the  differential  correction  procedure  will  be 
limited  to  an  improvement  of  the  6  Keplerian  elements.  This  mode  of  operation 
will  result  in  a  significant  saving  of  computer  time,  since  the  partial  derivatives 
of  present  position  with  respect  to  density  are  unnecessary. 

Parameters  33  and  34  al.'ow  the  user  to  specify  a  linear  function  for  the 
assumed  center  of  the  diurnal  atmospheric  bulge  as  follows: 

bulge  center  =  parameter  33  +  parameter  34  times 
where  5Q  =  declination  of  the  sun. 

Parameter  33  is  assumed  to  be  in  degrees  while  parameter  34  is  a  pure 
number. 

The  solar  position  is  only  recomputed  when  the  difference  between  the 
present  time  and  the  time  at  the  previous  computation  of  solar  position  is  more 
than  3  hours  as  indicated  by  parameter  35. 

Yhe  program  initially  attempts  to  process  epochs  which  span  9  hours  as 
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specified  by  parameter  40.  If  unsuccessful,  the  epoch  length  is  successively 
increased  by  an  amount  specified  by  parameter  41  until  the  program  succeeds 
in  processing  an  epoch  or  the  maximum  allowable  epoch  (as  specified  by  param¬ 
eter  43)  is  exceeded.  Note  that  an  unreasonably  small  value  for  parameter  40 
will  result  in  a  considerable  waste  of  computer  time.  Factors  to  be  considered 
in  choosing  an  appropriate  value  for  parameter  40  are  the  area-to-mass  ratio 
and  (perigee)  height  of  the  satellite,  the  observation  quality,  and  the  desired 
accuracy  of  the  density  result. 

The  mean  anomaly  residual  calculation  (which  is  used  by  the  preliminary 
adjustment  procedure)  may  at  times  give  slightly  inaccurate  results.  This 
deficiency  may  be  overcome  by  an  iterative  procedure.  Hence,  the  need  for 
parameter  47.  This  parameter  is  automatically  increased  by  1  for  the  first 
epoch  to  be  processed  and  for  any  epoch  which  follows  an  unsuccessful  dif¬ 
ferential  correction  attempt.  The  entire  preliminary  adjustment  procedure 
may  be  omitted  by  setting  parameter  47  to  -1. 

If  parameter  52  is  greater  than  zero,  the  results  of  each  iteration  of  the 
differential  correction  procedure  will  be  printed.  Otherwise,  a  printout  will 
occur  only  during  the  final  iteration. 

Parameter  55  is  automatically  decreased  by  1  for  any  epoch  in  which  the 
preliminary  adjustment  procedure  has  operated  successfully. 

The  differential  correction  procedure  is  considered  to  have  diverged, 
if  on  a  given  iteration  the  value  of  sigma  (standard  error  of  the  correction)  is 
five  times  the  sigma  of  the  previous  iteration  as  indicated  by  parameter  71. 

The  differential  correction  procedure  is  considered  to  have  converged  if 
2  successive  sigmas  agree  to  1%  as  indicated  by  parameter  72.  If  the  proce¬ 
dure  reaches  the  last  allowable  iteration  without  diverging  or  converging,  the 
procedure  is  considered  to  have  converged  if  the  last  2  sigmas  agree  to  10% 
as  indicated  by  parameter  73. 

Parameter  76  is  used  by  the  program  for  updating  right  ascension  and 
declination  measurements  from  a  given  equinox  to  a  common  equinox.  The 
program  initially  sets  parameter  76  equal  to  January  1  of  the  year  closest 
to  the  time  of  the  starting  elements  as  specified  on  the  starting  elements 
cards. 
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Parameter  77  is  set  by  the  program  to  a  height  (to  the  nearest 
ten  kilometers)  close  to  perigee  as  specified  on  the  starting  elements 
cards. 

Parameter  76  and  77  may  be  modified  by  the  user  in  the  usual 
manner. 

The  parameter  table  contains  all  zonal,  tes  serai,  and  sectorial 
harmonics  through  8,8.  Values  shown  in  Table  A  -  1  may  be.  modified 
by  means  of  the  following  two  formulas  which  give  the  correspondence 
between  a  parameter  number  and  a  given  harmonic. 

C™  is  placed  in  parameter  number  [80  +  +  m-2] 

S™  is  placed  in  parameter  number  [308  +  m  -  1] 

All  harmonic  coefficients  are  assumed  to  be  unnormalized. 

The  present  version  of  the  LGNDR  routine  which  evaluates  the 
associated  Legendre  functions  contains  all  functions  up  through  (8,8). 
Additional  functions  up  through  (20,  20)  are  available  and  may  be  easily  in¬ 
serted  into  the  LGNDR  deck.  Cm  or  Sm  requires  Pmand  pm+^  for  m  <  n 

n  n  n  n 

.  , 

and  P  for  m  =  n. 
n 

A.  6  DRIVE  Routine 

The  DRIVE  routine  calls  in  the  main  control  routine  ADMON,  which 
in  turn  activates  lower  level  routines  as  described  in  section  4.  The  logical 
numbers  of  the  system  input  and  output  tapes,  namely  5  and  6  respectively, 
are  transmitted  by  DRIVE  as  calling  sequence  arguments  to  ADMON.  An 
additional  degree  of  flexibility  has  been  incorporated  into  the  DRIVE  pro¬ 
gram  which  allows  the  user  to  generate  any  of  the  three  required  data  tapes 
at  run  time.  The  use  of  this  feature  is  described  in  paragraph  A.  7  of  this 
a  ppendix. 

A.  7  Data  Deck  Setup  for  Running  GADNIP 

Figure  A -l  shows  a  typical  data  deck  setup  for  three  CADNIP  runs. 
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In  this  figure,  capital  letters  are  used  for  cards  whose  type  is  determined 
by  a  specific  code  as  described  in  paragraph  A.  4  of  this  appendix.  The 
types  of  all  other  cards  appearing  in  this  figure  are  determined  by  their 
relative  positions  in  the  data  deck. 

The  first  card  shown  in  figure  A.  1  is  read  in  by  the  DRIVE  routine. 
This  card  contains  three  numbers  which  are  read  according  to  the  follow¬ 
ing  FORMAT  statement. 

FORMAT  (314) 

The  three  numbers  on  this  card  indicate  respectively  whether  the 
required  station,  density,  and  observation  tapes  are  available  from  a  previous 
operation,  or  are  to  be  generated  at  run  time.  Thus  in  the  example  of 
figure  A-l  ,  the  user  intends  to  generate  at  run  time  the  station  and  obser¬ 
vation  tapes,  while  a  density  tape  is  available  from  a  previous  operation. 

The  logical  numbers  of  tapes  generated  at  run  time  are  defined  by  DRIVE 
to  be  the  same  as  those  assigned  by  ADMON  in  the  parameter  table,  namely 
1,  2  and  3.  Thus,  in  this  example  25  cards  will  be  copied  onto  logical 
tape  I  and  the  next  50  cards  will  be  copied  onto  logical  tape  3.  Cards  used 
to  generate  data  tapes  in  this  fashion  must  be  identical  in  form  and  in  order 
to  cards  used  in  the  preparation  of  off-line  data  tapes.  Thus  the  25  station 
cards  must  include  all  cards  which  comprise  the  three  sections  of  a 
static. 'ii  tape  along  with  any  required  blank  cards  as  described  in  paragraph 
A.  2  of  this  appendix.  Similarly,  the  blank  card  which  follows  the  last  ob¬ 
servation  on  an  observation  tape  will  be  the  last  card  of  the  50  observation 
cards. 

After  generating  the  required  data  tapes,  the  DRIVE  routine  calls 
ADMON  which  in  turn  reads  and  processes  all  cards  up  to  and  including 
the  first  run  card.  When  the  time  interval  appearing  on  this  run  card 
has  been  fully  processed,  ADMON  reads  and  processes  the  next  set  of 
cards  up  to  and  including  the  second  run  card. 

Since  no  CHANGE  cards  (which  may  alter  the  logical  numbers  of 
data  tapes  to  be  processed)  appeared  in  either  of  the  first  two  CADNIP 
runs,  the  same  data  tapes  will  be  used  for  both  runs.  Thus  it  is  possible 
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to  use  run  time  generated  tapes  in  a  number  of  runs. 

The  END  card  which  is  then  read  by  ADMON  effectively  returns 
control  to  the  DRIVE  routine.  At  this  point,  DRIVE  reads  the  card 
containing  all  zeros  (indicating  that  no  data  tapes  are  to  be  generated) 
and  immediately  calls  ADMON.  ADMON  then  reads  and  processes  the 
third  CADNIP  run,  and  upon  encountering  the  last  END  card  again  re¬ 
turns  control  to  DRIVE.  The  entire  procedure  is  terminated  by  an  end-of- 
file  condition  encountered  by  DRIVE. 

It  will  be  observed  that  the  first  END  card  and  the  following  card 
containing  all  zeros  which  appear  in  figure  A-l  are  unnecessary  and  could 
have  been  omitted,  in  which  case  all  three  runs  would  have  been  processed 
in  succession  by  ADMON  without  returning  control  to  DRIVE  after  the 
second  run.  Note,  however,  that  it  would  be  improper  to  remove  only 
one  of  these  two  cards. 
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Introduction 

Satellite  positions  are  predicted  in  the  CADNIP  program  by  a  numerical 
integration  of  the  equations  of  motion  in  rectangular  coordinates,  expressed  in 
the  form  of  six  simultaneous  first-order  equations.  The  integration  algorithm 
is  a  predictor-corrector  method,  using  a  sixth-order  Adams-Bashforth  pre¬ 
dictor,  a  sixth  -  order  Adams  -  Moulton  corrector  (which  is  applied 
only  once),  and  modifiers  of  so-called  Hamming  type  to  improve  both  predictor 
and  corrector.  Starting  values  are  generated  by  the  classical  fourth-order 
Runge-Kutta  method. 

Stability,  accuracy,  and  simplicity  were  the  principal  criteria  conside  red 
in  choosing  the  integration  method  to  be  implemented  in  CADNIP.  For  a  dis¬ 
cussion  of  the  relative  merits  of  the  various  approac'  s  to  the  numerical  in¬ 
tegration  of  satellite  trajectories,  the  reader  is  referred  to  Conte. 

For  the  sake  of  notational  uniformity,  the  equations  of  motion  will  be 


expressed  in  vector  form  as 


Y  -  F(Y,t) 


and  subscripts  will  be  used  to  designate  values  at  particular  instants  of  time, 


Y.  =  Y(t.)  and  F.  =  F(Y.,t.) 

11  i  11 

1.  The  Runge-Kutta  Starter 

The  multistep  predictor-corrector  method  used  in  CADNIP  requires  six 
previous  values  of  the  right-hand  sides,  F,  and  of  the  state  vector,  Y,  in  or¬ 
der  to  generate  one  new  step  of  the  solution.  At  the  beginning  of  an  integration, 
these  starting  values  are  generated  by  a  self-starting  single  -  step  method,  the 

12 

classical  fourth-order  Runge-Kutta  algorithm.  (See,  for  example,  Hildebrand  , 
p.  237;  Henrici  ,  pp.  121-122.) 

There  is  a  provision  in  CADNIP  for  the  use  of  three  different  step  sizes, 
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depending  on  the  purpose  for  which  the  numerical  integration  is  being  done: 

.  ephemeris  generation 

.  evaluation  of  partial  derivatives  with  respect 
to  initial  state 

.  evaluation  of  partial  derivatives  with  respect 
to  density 

All  three  of  these  parameters  are  currently  set  at  15  seconds.  Similarly, 
the  ratio  between  the  Runge-Kutta  step  size  and  the  step  size  to  be  used  by 
thepredictor  -corrector  algorithm  may  also  be  chosen  individually  for  each  of 
the  three  cases  listed  above;at  present,  however,  each  of  these  three  ratios 
is  set  at  4.  Thus  a  total  of  20  Runge-Kutta  steps  are  taken  in  order  to  gen¬ 
erate  the  required  set  of  starting  values. 

The  reader  may  consult  section  A.  5  for  information  on  changing  these 
parameters. 

2.  The  Predictor-Corrector  Method 

The  basic  algorithm  consists  of  four  steps: 

.  Predict 
.  Modify  predictor 
.  Correct 
.  Modify  corrector 

at  the  end  of  which  a  new  state  vector  has  been  computed.  The  basic  for¬ 
mulas  are  derived  in  sections  2.  1  and  2.  2,  below;  a  concise  summary  and 
flow  chart  are  presented  in  section  2.  3. 

2.  1  Basic  Predictor  -  Cor  rector  Formulas 

For  the  sake  of  simplicity,  the  formulas  are  derived  here  in  scalar  form. 
The  extension  to  vector  variables  is  immediate.  References  12,  13  and  14  may 
be  consulted;  however,  only  reference  14  gives  the  explicit  form  of  formula 
(B.  5)  correctly. 

2,  1,  1  The  Predictor 

We  write  y  =  y  +  Ay,  where  Ay  represents  the  change  in  y  between  the 

n+ 1  n 

times  t  =  t  and  t  =  t  That  is, 

n  n+  1 
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Ay  =  I  n+1  (iyj  dt 
dt 


Since  the  differential  equation  is  of  the  form 


±L  _ 

dt 


it  follows  that 


f(y.t) 


rVn 


^  y  =  i  f(y,  t)  dt  where  t  =  t  +  h 

7  /  n+ 1  n 


(B.  2) 


Of  the  many  forms  for  representing  a  collocation  polynomial  approximating 
f ( y ,  t),  we  choose  one  employing  backward  differences: 

V7  f  =  f  -  f  . 

V  n  n  n  ••  1 

V2^  =V(Vf  )  =Vf  H7E  ,  »  £  -f  ,  -  (f  ,-f  ,) 

n  n  n  n-1  n  n-1  n-1  n-2 


=  f  -2f  +  f 

n  n-1  n-2 


etc . 


In  terms  of  these  backward  differences,  the  polynomial  P(t),  agreeing  with 

f(v,t)  at  the  points  t,t  ,  ,  t  ,, 

n  n-1  r.  -  k 

can  be  written  in  the  form: 


P(t)  =  f  +  sVf  +  6(8+1)  r/f  +  .  .  .  +  s(s+l)..  .(s  +  k-1)  K 
n  n  2!  n  k!  n 


(B,  3) 


.  t-t  ,  . 

where  s  =  n  or  t  =  t  +  sh 

IT 


Then 


n+1 


1 

P(t)  dt  =  I  J^fn+  s\7fn+-  •  .J  (hds) 


fl*'n+  1 


t  P(t)dt  =  h^  cw7“fn 
oc=  0 


B  -  3 


where 


ci  i0  =  i 

s-^-  =  rf 

T^d-| 

C  r1  818+1118  +  2118  +  3)  _  2FJ. 

4  ^  4!  as  720 

_  f*  s(s  +  l)(8  +  2)(s  +  4)  _  95 

5  JQ  51  "  288  ~ 


The  final  form  is 


etc. 


Vi  ■  yn  *  h«„  4vf„ +  if  v\  *  I  v\  *  Mi  \\  *  2§  v\> 


+  Remainder  term 


(B.4) 


if  fifth  differences  are  retained.  This  is  the  Adams  -  Bashforth  formula 
in  the  backward-difference  form. 

To  transform  from  backward  difference  to  ordinate  form,  we  use  the 
relations 


Vkfn  =  V(l)fn-l+(2)fn-2 


-  ....  +  (-1)  f 


n-k 


Considering,  in  particular,  formula  (B.4) 


y  , ,  =  i  y  +h(f  +  . . .  +  ||-  y5f  ) 

n+ 1  ^  n  n  288  v  n 


+  Remainder 


we  obtain,  upon  carrying  out  the  substitutions  and  collecting  terms, 
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Y  ,  ■  y  +  7777  (4277  f  -  7923f  +  9982f  -  7298£  ,  + 

7n+l  7n  1440  n  n-1  n-2  n-3 


475  £  .)  +  Remainder 

n- 5 


This  is  the  predictor  formula  used  in  CADNIP. 

2.1,2  The  Corrector 

Another  way  of  defining  a  collocation  polynomial  P(t)  to  ap; 

f(y.t)  would  be  to  require  agreement  at  the  points  t  , ,  t  ,  ... 

n  I  i  n 

P(t)  would  then  be  of  the  form 

* 

*  *  s  (  R  J  1  'I  2 

p(*)afn+l  +  i^£n+l  it  Vfn+1  +  - 


>;<  t  - 1  * 

where  s  =  n+1  or  t  =  t  ,  +  s  h 

— : -  n+ 1 


Now 


,t 
"  t 


rn+ip*(.)d« 


n 


•  n 


,ntl  +  *V'n+ !+"•] 


n+ 1  v"  ** 

,  p  <'»<“=  hZ  ^Vfn+1 


where 


C  =  1 


*  *  1 

s  ds  =  -  — 


-i: 


a0  *  *  , 

r  s  (s  4i) 

j.  1  2! 


■ds  =  -  — 


1_ 

12 


■  r 


s  (s  +  1)  (5  +  2) 
3! 


ds  =  -  — 


1_ 

24 


etc. 


Retaining  fifth  differences,  as  before,  we  obtain  the  formula 


2877f  , 

n-4 

(B.  5) 

:roxir.iate 
’  ‘n-k+l* 

(B.6) 

(h  ds  ) 


3  -  5 


yn+l 


n+1  2  v  n+1  12  v  n+1  24  v  n+1  720  v  n+1  160 v 


- 


+  Remainder  term 
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which  is  ihe  Adams -Moulton  formula  in  backward-difference  form.  Since 

(B .  7)  expresses  y  ,  in  terms  of  f  , ,  which  is  itself  a  function  of  y  , 

'  n+1  n+1  n+l 

the  formula  is  implicit.  It  must  be  used  in  conjunction  with  a  predictor  of 
explicit  type  (e.g.  formula  B.4)  which  supplies  a  value  of  f 

The  transformation  to  ordinate  form  is  carried  out  as  before,  with 


the  result 


y  ,  =  y  +  T77o(475  f  +  1427  f  -  798f  .  +482f  _  -  173  f  +  27f 

7nt  1  n  1440  n+1  n  n- 1  n-2  n-3 


+  Remainder  term  (B.8 

This  is  the  corrector  formula  used  in  CADN1P. 


2.  1.3  The  Remainder  Terms 

The  local  quadrature  error  introduced  by  using  a  finite  number  of  dif¬ 
ferences  has  as  its  principal  part  the  first  neglected  term.  Thus,  f<~  the 
predictor  formula. 


y 


n+  1 


i  y  +  h(f  +  .  .  .  + 
n  n 


95 

288 


t  h( 


19087 

60480 


)  V6f 


as  h - ►  0 


For  any  n,  the  backward  difference  may  be  expressed  in  the  form 


Thus 


V  f  --  h  2_ 

dtn 

hy*f  =  h7.d_  f<t  ) 
v  n  6  n' 

dt 


h"dn 


Since 


6 

h  V  f  = 
n 


h 


7  [7] 

y 


Therefore  the  principal  error  term  for  the  predictor  is 

o  ■  .  19087  7  £71 

Remainder  =  +  h  y 
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The  principal  error  term  for  the  corrector  is  derived  in  exactly  the  same 


863  .  7  [7] 

Remainder  =  -  -ggjgg  ”  h  y 


(B.  io) 


2.  2  The  Modifiers 


If  y.  .  designates  the  true  value  of  y(t.  , )  then,  in  the  limit  as  h  — *■  0, 
l+l  l+l 


(p)  19087  ,  7  C7) 

yi+l  =  Vi+1  +  60480  h  y 


(B.  1  1) 


(pO 

where  y.  +  j  is  the  right-hand  side  of  (B.  5)  excluding  the  remainder 
term.  Similarly,  from  (B.  8)  and  (B.  10) 


(c)  863  .  7  17] 

Vi  =  vTi  '63«5‘hy 


(B.  12) 


Subtracting  (B.  12)  from  (B.  11),  we  obtain  in  the  limit  (dropping  the 
subsc ript  i+ 1) 


=  r<P>_  (C)1+  19087  +J63h7  17]  = 

L y  y  J  60480  v 


(p)  (c)  + 

y  -  y 


’9930  ,7  [7 ] 
60480  b  y 


,  7  [7]  60480  (c)  (P)  ' 

h  r  =  y  -  y  j 


This  approximation  for  h*y  L  j  permits  us  to  estimate  the  error  terms  (B.  9) 
and  (B.  10).  Applying  these  estimates  as  corrections,  or  modifications,  to 
the  predicted  and  corrected  values,  we  obtain,  finally,  thr.  modified  predictor 
and  corrector  values: 


(p.m)  (p)  19087  7  i.7j  (p)  19087  .  60480  ;  (c)  <p)  ■ 

y  =  y  +  i  h  y  =  y  T  mrr.  at:  y  -  y 


60480 


60-130  19Q50 


fp.m)  r  (P)  +  i9087  r  (c)  _  (pfj 
f  y  19950  Ly  y  J 


previous  step 


y(p)  +  (  0.956741854  )  j  y(c)  -  y(p)]  . 

L  J  previous  step 


(c,m)  (c)  863  .  60480  !  (c)  (p)!  (c)  863  \  (c)  <p) ! 

y  "  y  ‘60480  19950  [y  y  j  =  y  '  19950  iy  y  J 

y(c,m)  _  [‘19087  y(c>  +  863  y(p)  |  ^19950... 

1  5 

The  use  of  modifiers  of  this  type  is  often  attributed  to  Hamming  , 
although  the  basic  idea  goes  back  at  least  to  Richardson^. 


2.  3  Flow  Chart 

Assume  that  six  previous  state  vectors  Y  and  right  hand  sides  F  are 
available,  either  from  previous  predictor -corrector  steps  or  (at  the  begin¬ 
ning  of  an  integration)  from  the  Runge-Kutta  starting  algorithm. 


Time 

s 


State  vector 


Right-hand  sides 
F1 


t . 


Alsu  required  are  the  intermediate  predictor  and  corrector  values  at  the 
sixth  step,  YP,  and  YC  .  Then  the  following  flow  chart  describes  the 

D  6 

foul  operations  leading  to  the  values 
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Time 


State  vector 


Rieht-hand  sides 


rr'T  r'l! 


Notation 


YP 

YPM 

YC 

Y 

FP 

F 

h 


predicted  values 
modified  predicted  values 
corrected  values 

modified  corrected  values  (final  values) 
right-hand  sides  computed  using  YPM 
right-hand  sides  computed  using  Y  (final  values) 
step  size 


Constants 


i  =  0 


2  3  4  5 


4277  -7923  9982  -7298  2877  -475 


(7  =  0.  956741854 


i  =  0  1  2  3  4  5 

=  475  1427  -798  482  -173  27 

i 

A  =  19087 

/v  =  863 

V  =  19950 


1 
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3.  Testing  of  the  Predictor-Corrector  Algorithm 


Assurance  that  the  numerical  integration  algorithm  operates  at 
the  required  level  of  accuracy  is  guaranteed  by  a  series  of  tests  performed 
both  before  and  after  the  algorithm  was  incorporated  into  CADNIP. 

1.  To  verify  the  accuracy  of  the  basic  predictor-corrector 
formulas,  tests  were  run  using  a  Runge-Kutta  formula,  a  fourth-order 
Adams -Bashforth-Moulton  procedure,  and  the  presently  implemented 
sixth-order  Adams-Bashforth-Moulton  method,  to  generate  positions  in 
Keplerian  orbits  of  various  eccentricities, 

2.  The  integration  method  with  modifiers  was  compared 
with  the  basic  unmodified  method,  demonstrating  the  increased  accuracy 
resulting  from  the  use  of  modifiers.  No  loss  of  stability  was  observed, 

3.  The  truncation  error  growth  as  a  function  of  step  size 
was  found  to  conform  with  theory.  Numerical  evidence  indicated  that  the 
growth  of  round  off  error  was  insignificant, 

4.  Runs  were  made,  using  typical  satellite  orbits,  to  de¬ 
termine  empirically  the  most  appropriate  step  sizes  for  the  Runge-Kutta 
and  predictor-corrector  algorithms. 

5.  With  the  predictor-corrector  method  incorporated  op¬ 
erationally  into  the  CADNIP  program,  test  runs  were  made  using  actual 
satellite  data: 

i)  using  Runge-Kutta  alone 

ii)  using  Runge-Kutta  only  as  a  starter  with  the 
predictor-corrector  method  being  used  to  gen¬ 
erate  the  orbit. 

The  results  agreed  with  one  another,  and  both  reproduced  the  satellite  mo¬ 
tion  to  within  an  acceptable  level  of  residuals. 

4.  Stability  of  the  Method 

As  a  practical  matter,  the  staoility  of  the  predictor-modifier- 
corrector -modifier  algorithm  has  been  established  empirically  by  the  tests 
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described  in  the  previous  section,  as  well  as  by  the  absence  of  symptoms 
of  instability  during  the  extensive  production  runs  using  the  CADNIP  prog¬ 
ram. 

From  a  theoretical  point  of  view,  the  stability  of  a  class  of  algo¬ 
rithms  including  the  method  used  in  CADNIP  has  been  investigated  by 
Abdel  Karim, ^  whose  principal  result  is  a  theorem  giving  necessary  and 
sufficient  conditions  for  stability  in  terms  of  the  positive -definitenes s  of 
certain  Hermitian  forms  associated  with  the  system  of  differential  equations. 
Unlike  earlier  asymptotic  stability  theorems  which  apply  only  in  the  limit 
as  h  —*  0,  the  criteria  of  Abdel  Karim  are  functionally  dependent  on  h  and 
can  be  used,  therefore,  to  define  stability  boundaries.  It  is  clear,  however, 
fur  the  particular  method  being  considered  here,  that  values  of  h  large  e- 
nough  to  induce  instability  do  not  arise  in  the  CADNIP  application,  where 
the  step  rize  is  limited  by  the  relatively  stringent  propagated  truncation 
error  requirement. 
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